home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_6 / a6_2.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.5 KB  |  117 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 6.2 (Differentiation Using Extrapolation).
  9. % Section    6.1,    Approximating the Derivative, Page 327
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program finds a numerical approximation for f'(x0).
  13.  
  14. % The method employed is Richardson`s extrapolation.
  15.  
  16. %    The function f(x) is stored in the M-file  f.m.
  17. % function z = f(x)
  18. % z = cos(x);
  19.  
  20. delete f.m
  21. diary f.m; disp('function z = f(x)');...
  22.            disp('z = cos(x);');...
  23. diary off;
  24.  
  25. f(0); % Remark. f.m and diffext.m are used for Algorithm 6.2
  26.  
  27. pause % Press any key to continue.
  28.  
  29. clc;
  30. %    Place abscissa for differentiation in x0.
  31.  
  32. % Place the endpoints interval [a,b] about x0 in a and b.
  33.  
  34. x0 = 0.8;
  35.  
  36. a  = 0;
  37.  
  38. b  = pi/2;
  39.  
  40. y0 = f(x0);
  41.  
  42. pause    % Press any key to graph the function.
  43.  
  44. clc; clg;
  45. a = 0;
  46. b = 1.6;
  47. c = 0;
  48. d = 1.1;
  49. h = (b-a)/100;
  50. Xs = a:h:b;
  51. Ys = f(Xs);
  52. axis([a b c d]);...
  53. plot(x0,y0,'or',Xs,Ys,'g');...
  54. hold on;...
  55. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  56. xlabel('x');...
  57. ylabel('y');...
  58. title('y = f(x) and the given point');...
  59. grid;...
  60. axis;...
  61. hold off;...
  62. shg; pause    % Press any key to continue.
  63.  
  64. pause    % Press any key to find the derivative.
  65.  
  66. clc; clg;
  67.  
  68. %    Place abscissa for differentiation in x0.
  69.  
  70. % Place the tolerance for the error in  delta.
  71.  
  72. % Place the tolerance for the relative error in  toler.
  73.  
  74. x0 = 0.8;
  75.  
  76. delta = 1e-10;
  77.  
  78. toler = 1e-10;
  79.  
  80. [D,err,relerr,n] = diffext('f',x0,delta,toler);
  81.  
  82. pause % Press any key to see the solution.
  83.  
  84.  
  85.  
  86. X1 = [a b];
  87. C = [D(n,n) f(x0)];
  88. Y1 = polyval(C,X1-x0);
  89. axis([a b c d]);...
  90. plot(x0,y0,'or',Xs,Ys,'-g',X1,Y1,'-r');...
  91. hold on;...
  92. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  93. xlabel('x');...
  94. ylabel('y');...
  95. title('The tangent line has slope m = f`(x0).');...
  96. grid;...
  97. axis;...
  98. hold off;...
  99. shg; pause    % Press any key to continue.
  100.  
  101. format long;
  102. Mx0 = 'Approximating the derivative with Richardson`s extrapolation.';
  103. Mx1 = 'Difference table of approximate derivatives for ';
  104. Mx2 =  ['f`(',num2str(x0),')'];
  105. Mx3 = [Mx1,Mx2];
  106. Mx4 = 'The approximate value of ';
  107. Mx5 = [Mx4,Mx2,' = '];
  108. Mx6 = 'An estimate for the error bound is:';
  109. Mx7 = ' approx. derivative   +- error bound';
  110. clc,echo off,diary output,...
  111. disp(''),disp(Mx0),disp(''),disp(Mx3),disp(D),...
  112. disp(''),disp(Mx5),disp(D(n,n)),...
  113. disp(Mx6),disp(Mx7),disp([D(n,n) err]),diary off,echo on
  114.  
  115.  
  116.  
  117.